package scu.maqiang.numeric;

import scu.maqiang.csparsej.tdouble.Dcs_common;
import scu.maqiang.csparsej.tdouble.Dcs_entry;
import scu.maqiang.csparsej.tdouble.Dcs_util;
import scu.maqiang.cxsparsej.DZcs_common;
import scu.maqiang.cxsparsej.DZcs_entry;
import scu.maqiang.cxsparsej.DZcs_util;

import java.util.*;
import java.util.stream.IntStream;

import static java.lang.Math.*;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;

/**
 * 复稀疏矩阵类
 * 
 * @author 马强
 * @version 1.0
 * @since 1.0
 *
 */
public class SCMatrix {

	private int m = 0;
	private int n = 0;
	public int nnz = 0;
	private boolean sorted = false;
	public ArrayList<SCEntry>[] mat;

	/**
	 * 构造函数，指定稀疏矩阵行数与列数产生空的稀疏矩阵
	 * 
	 * @param m 矩阵行数
	 * @param n 矩阵列数
	 */
	@SuppressWarnings("unchecked")
	public SCMatrix(int m, int n) {
		this.m = m;
		this.n = n;
		this.nnz = 0;
		mat = new ArrayList[m];
		for (int i = 0; i < m; i++) {
			mat[i] = new ArrayList<SCEntry>();
		}
	}

	/**
	 * 构造方型复稀疏矩阵
	 * 
	 * @param n 矩阵阶数
	 */
	public SCMatrix(int n) {
		this(n, n);
	}

	/**
	 * 通过给定三元组构造复稀疏矩阵
	 * 
	 * @param i     矩阵行指标
	 * @param j     矩阵列指标
	 * @param value 矩阵元素值
	 * @param row   矩阵行数
	 * @param col   矩阵列数
	 */
	public SCMatrix(int[] i, int[] j, Complex[] value, int row, int col) {
		this(row, col);
		for (int ii = 0, nn = i.length; ii < nn; ii++) {
			if ((value[ii].norm()) == 0.0) {
				mat[i[ii]].add(new SCEntry(j[ii], value[ii]));
				this.nnz++;
			}
		}
	}

	/**
	 * 通过给定三元组构造复稀疏矩阵
	 * 
	 * @param i     矩阵行指标
	 * @param j     矩阵列指标
	 * @param value 矩阵元素值ֵ
	 * @param size  矩阵阶数
	 */
	public SCMatrix(int[] i, int[] j, Complex[] value, int size) {
		this(i, j, value, size, size);
	}

	/**
	 * 复制构造函数
	 * 
	 * @param sm 复稀疏矩阵
	 */
	@SuppressWarnings("unchecked")
	public SCMatrix(SCMatrix sm) {
		this.m = sm.m;
		this.n = sm.n;
		mat = new ArrayList[m];
		for (int i = 0; i < this.m; i++) {
			mat[i] = new ArrayList<SCEntry>();
			for (SCEntry ele : sm.mat[i]) {
				mat[i].add(new SCEntry(ele));
				this.nnz++;
			}
		}
	}

	/**
	 * 分别传入两个实矩阵作为复矩阵的实部与虚部
	 * @param mat_real
	 * @param mat_imag
	 */
	public SCMatrix(SRMatrix mat_real, SRMatrix mat_imag) {
		this.m = mat_real.getM();
		this.n = mat_real.getN();
		assert(m == mat_imag.getM());
		assert(n == mat_imag.getN());
		mat = new ArrayList[m];
		for(int i = 0; i < m; i++) {
			mat[i] = new ArrayList<SCEntry>();
			for(SREntry ele: mat_real.mat[i]) {
				setElementFast(i, ele.col, new Complex(ele.val, 0.0));
			}
			for(SREntry ele: mat_imag.mat[i]) {
				addElement(i, ele.col, new Complex(0.0, ele.val));
			}
		}

	}
	
	@SuppressWarnings("unchecked")
	public SCMatrix(SRMatrix sr) {
		this.m = sr.getM();
		this.n = sr.getN();
		mat = new ArrayList[m];
		for (int i = 0; i < this.m; i++) {
			mat[i] = new ArrayList<SCEntry>();
			for (SREntry ele : sr.getRow(i)) {
				mat[i].add(new SCEntry(ele.col, new Complex(ele.val, 0.0)));
				this.nnz++;
			}
		}		
	}

	/**
	 * ���ⲿ�ļ��ж�ȡϡ�����
	 * 
	 * @param matFile �ⲿ�ļ���
	 */
	@SuppressWarnings("unchecked")
	public SCMatrix(String matFile) {
		Scanner sc = null;
		try {
			sc = new Scanner(new File(matFile));
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		}

		int m = sc.nextInt();
		int n = sc.nextInt();
		int nnz = sc.nextInt();
		String isSorted = sc.next();
		sc.next();
		this.m = m;
		this.n = n;
		mat = new ArrayList[m];
		for (int i = 0; i < this.m; i++) {
			mat[i] = new ArrayList<SCEntry>();
		}
		for (int i = 0; i < nnz; i++) {
			int iIdx = sc.nextInt();
			int jIdx = sc.nextInt();
			double realValue = sc.nextDouble();
			double imageValue = sc.nextDouble();
			mat[iIdx - 1].add(new SCEntry(jIdx - 1, realValue, imageValue));
		}
		if (isSorted.trim() == "T") {
			this.sorted = true;
		} else {
			this.sorted = false;
		}
	}

	/**
	 * 将矩阵保存到一个文件中
	 * 
	 * @param matFile 文件函数名
	 */
	public void save(String matFile) {
		sort();
		boolean isHermitian = isHermitian();
		String isSym = isHermitian ? "T" : "F";
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(matFile))) {
			bw.write(m + "\t" + n + "\t" + nnz + "\tT\t" + isSym + "\n");
			for (int i = 0; i < m; i++) {
				for (SCEntry ele : mat[i]) {
					bw.write((i + 1) + "\t" + (ele.col + 1) + "\t" + ele.val.real + "\t" + ele.val.image + "\n");
				}
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	/**
	 * 清除矩阵中所有元素
	 */
	public void clear() {
		for (int i = 0; i < m; i++) {
			mat[i].clear();
		}
		nnz = 0;
		sorted = false;
	}

	public void setDiag(Complex val) {
		for (int i = 0; i < m; i++) {
			mat[i].clear();
			mat[i].add(new SCEntry(i, val));
		}
		nnz = m;
		sorted = true;
	}

	public void exchangeRow(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		ArrayList<SCEntry> temp = new ArrayList<>();
		for (SCEntry ele : mat[i]) {
			temp.add(new SCEntry(ele.col, ele.val));
		}
		mat[i].clear();
		for (SCEntry ele : mat[j]) {
			mat[i].add(new SCEntry(ele.col, ele.val));
		}
		mat[j].clear();
		for (SCEntry ele : temp) {
			mat[j].add(new SCEntry(ele.col, ele.val));
		}
	}

	public void exchangeCol(int i, int j, int nBand) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		ArrayList<SCEntry> iColTemp = new ArrayList<>();
		ArrayList<SCEntry> jColTemp = new ArrayList<>();
		int iIdxBegin = Math.max(0, i - nBand - 1);
		int iIdxEnd = Math.min(m, i + nBand + 1);
		for (int iIdx = iIdxBegin; iIdx < iIdxEnd; iIdx++) {
			for (int jj = 0, n = mat[iIdx].size(); jj < n; jj++) {
				SCEntry temp = mat[iIdx].get(jj);
				if (temp.col == j) {
					iColTemp.add(new SCEntry(iIdx, temp.val));
					mat[iIdx].remove(jj);
					break;
				}
			}
		}

		// System.out.println(MVO.toString(this.toDenseMatrix()));
		int jIdxBegin = Math.max(0, j - nBand - 1);
		int jIdxEnd = Math.min(m, j + nBand + 1);
		for (int jIdx = jIdxBegin; jIdx < jIdxEnd; jIdx++) {
			for (int ii = 0, n = mat[jIdx].size(); ii < n; ii++) {
				SCEntry temp = mat[jIdx].get(ii);
				if (temp.col == i) {
					jColTemp.add(new SCEntry(jIdx, temp.val));
					mat[jIdx].remove(ii);
					break;
				}
			}
		}

		// System.out.println(MVO.toString(this.toDenseMatrix()));

		for (SCEntry ele : iColTemp) {
			mat[ele.col].add(new SCEntry(i, ele.val));
		}

		for (SCEntry ele : jColTemp) {
			mat[ele.col].add(new SCEntry(j, ele.val));
		}
	}

	public Pair<Integer, Complex> maxNormColValDown(int i, int nBand) {
		int end = Math.min(i + nBand + 3, m);
		Complex maxval = getElement(i, i);
		int maxIdx = i;
		for (int ii = i + 1; ii < end; ii++) {
			Complex temp = getElement(ii, i);
			if (temp.norm() > maxval.norm()) {
				maxval = temp;
				maxIdx = ii;
			}
		}
		return new Pair<Integer, Complex>(maxIdx, new Complex(maxval));
	}

	/**
	 * 将稀疏矩阵转化为一个稠密矩阵
	 * 
	 * @return 二阶数组
	 */
	public Complex[][] toDenseMatrix() {
		Complex[][] result = new Complex[m][n];
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				result[i][ele.col].real = ele.val.real;
				result[i][ele.col].image = ele.val.image;
			}
		}
		return result;
	}

	/**
	 * �õ���i������ָ��Ϊj��SCEntry����
	 * 
	 * @param i �����i��
	 * @param j �����j��
	 * @return ����SCEntry����, ������col����Ϊj, valΪԪ�ص�ֵ,
	 *         ���ϡ������в����ڵ�i�е�j��Ԫ��, �򷵻�null
	 */
	public SCEntry getSCEntry(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		for (SCEntry ele : mat[i]) {
			if (ele.col == j) {
				return ele;
			}
		}
		return null;
	}

	/**
	 * �õ�����Ԫ�ص�ֵ
	 * 
	 * @param i Ԫ����ָ��
	 * @param j Ԫ����ָ��
	 * @return ������ڣ��򷵻�Ԫ�ص�ֵ�����򷵻�0
	 */
	public Complex getElement(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		Optional<SCEntry> element = mat[i].stream().filter(ele -> ele.col == j).findFirst();
		return element.isPresent() ? new Complex(element.get().val) : new Complex(0.0, 0.0);
	}

	/**
	 * 设置矩阵元素ֵ
	 * 
	 * @param i     元素行数
	 * @param j     元素列数
	 * @param value 元素值
	 * @return 返回对象本身，便于继续设置矩阵元素
	 */
	public SCMatrix setElement(int i, int j, Complex value) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		int size = mat[i].size();
		if (size == 0 && value.norm() > Constants.Er) {
			mat[i].add(new SCEntry(j, value));
			this.nnz++;
			return this;
		}

		for (int idx = 0; idx < size; idx++) {
			SCEntry se = mat[i].get(idx);
			if (se.col == j) {
				if (value.norm() > Constants.Er)
					se.val = value;
				else {
					mat[i].remove(idx);
					--(this.nnz);
				}
				return this;
			}
		}

		mat[i].add(new SCEntry(j, value));
		++(this.nnz);
		return this;
	}

	/**
	 * 快速设置矩阵元素
	 * 
	 * @param i     元素行数
	 * @param j     元素列数
	 * @param value 元素值
	 * @return 返回对象本身，便于继续设置矩阵元素
	 */
	public SCMatrix setElementFast(int i, int j, Complex value) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		mat[i].add(new SCEntry(j, value));
		++(this.nnz);
		return this;
	}

	/**
	 * ��ϡ�������Ӧ����ָ���ϼ�����ӦԪ��ֵ, ���ս�������������: 1. ���Ԫ��ֵΪ0�������κβ���
	 * 2. ���Ԫ��ֵ��0�����Ҿ���������ָ���д���Ԫ�أ���ԭ��Ԫ��ֵ��֮���, ���Ӻ����Ϊ0,
	 * ��ɾ����Ԫ�أ������0�����޸ĸ�Ԫ��ֵ 3.
	 * ���Ԫ��ֵ��0���������в����ڸ�Ԫ�أ���������������Ӹ�Ԫ��ֵ
	 * 
	 * @param i   Ԫ����ָ��
	 * @param j   Ԫ����ָ��
	 * @param val Ԫ��ֵ
	 * 
	 */
	public void addElement(int i, int j, Complex val) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		if (val.norm() > Constants.Er) {
			int size = mat[i].size();
			if (size == 0) {
				mat[i].add(new SCEntry(j, val));
				this.nnz++;
				return;
			}
			for (int idx = 0; idx < size; idx++) {
				SCEntry se = mat[i].get(idx);
				if (se.col == j) {
					double tempReal = se.val.real + val.real;
					double tempImage = se.val.image + val.image;
					if (sqrt(tempReal * tempReal + tempImage * tempImage) == 0.0) {
						mat[i].remove(idx);
						--(this.nnz);
					} else {
						se.val = new Complex(tempReal, tempImage);
					}
					return;
				}
			}
			mat[i].add(new SCEntry(j, val));
			++(this.nnz);
			return;
		}
	}

	public void addElement(int i, int j, double value) {
		this.addElement(i, j, new Complex(value, 0.0));
	}

	/**
	 * �ھ�����ɾ����i�е�i������Ԫ��ֵ���������ϡ����������з�0Ԫ�ص�����ָ���ǶԳƵ�
	 * 
	 * @param i ��ָ��
	 */
	public void removeSymmetricIJ(int i) {
		mat[i].forEach(ele -> {
			int col = ele.col;
			if (col != i) {
				removeElement(col, i);
			}
		});
		removeRow(i);
	}

	/**
	 * ȥ��ϡ�������ĳ����ָ���µ�Ԫ��
	 * 
	 * @param i Ԫ����ָ��
	 * @param j Ԫ����ָ��
	 * @return ���ɾ����Ԫ�ش��ڣ��򷵻�true; ���򷵻�false
	 */
	public boolean removeElement(int i, int j) {
		ParamCheck.checkBounds(i, 0, m);
		ParamCheck.checkBounds(j, 0, n);
		int size = mat[i].size();
		for (int idx = 0; idx < size; idx++) {
			if (mat[i].get(idx).col == j) {
				mat[i].remove(idx);
				--(this.nnz);
				return true;
			}
		}
		return false;
	}

	/**
	 * ɾ��ϡ�����ĳ������Ԫ��
	 * 
	 * @param i ������ָ��
	 */
	public void removeRow(int i) {
		ParamCheck.checkBounds(i, 0, m);
		this.nnz -= mat[i].size();
		mat[i].clear();
	}

	/**
	 * �õ���������
	 * 
	 * @return ��������
	 */
	public int getM() {
		return m;
	}

	/**
	 * �õ���������
	 * 
	 * @return ��������
	 */
	public int getN() {
		return n;
	}

	/**
	 * �õ������0Ԫ�ظ���
	 * 
	 * @return �����0Ԫ�ظ���
	 */
	public int getNNZ() {
		return nnz;
	}

	public int getNNZL() {
		int sum = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			sum += mat[i].stream().filter(ele -> ele.col < ii).count();
		}
		return sum;
	}

	public int getNNZU() {
		int sum = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			sum += mat[i].stream().filter(ele -> ele.col > ii).count();
		}
		return sum;
	}

	public Pair<Integer, Integer> getNNZLU() {
		int numL = 0;
		int numU = 0;
		for (int i = 0; i < m; i++) {
			int ii = i;
			numL += mat[i].stream().filter(ele -> ele.col < ii).count();
			numU += mat[i].stream().filter(ele -> ele.col < ii).count();
		}
		return new Pair<Integer, Integer>(numL, numU);
	}

	/**
	 * ʵ��ϡ��������Ԫ��, �����ǰ����ΪA, �������Ϊalpha��B, ���յ÷��ؾ���ֵCΪ C �� A +
	 * alpha * B.
	 * 
	 * @param alpha ʵ��
	 * @param       B����Ӿ���
	 * @return�����ؾ���ֵΪC = A + alpha * B
	 */
	public SCMatrix add(double alpha, SCMatrix B) {
		SCMatrix result = new SCMatrix(this);
		int M = result.m;
		boolean flag = false;
		if (abs(alpha) < Constants.Er) {
			return result;
		}

		for (int i = 0; i < M; i++) {
			int mat2Size = B.mat[i].size();
			for (int j = 0; j < mat2Size; j++) {
				// flag��¼�Ƿ��ڵڶ��������д��ڸ�ָ��
				flag = false;
				SCEntry e2 = B.mat[i].get(j);
				int resultSize = result.mat[i].size();
				for (int k = 0; k < resultSize; k++) {
					SCEntry e1 = result.mat[i].get(k);
					// Ѱ���������е���ͬ�����Ƿ�����ͬ����ָ��
					if (e1.col == e2.col) {
						double temp2Real = e1.val.real + alpha * e2.val.real;
						double temp2Image = e1.val.image + alpha * e2.val.image;
						// �������ͬ��ָ�꣬����֮������ֲ�Ϊ0,���޸���ָ����ֵ
						if (sqrt(temp2Real * temp2Real + temp2Image * temp2Image) > Constants.Er) {
							e1.val = new Complex(temp2Real, temp2Image);
						} else {
							// ������Ϊ0,������Ԫ��
							result.mat[i].remove(k);
							result.nnz--;
						}
						flag = true;
						break;
					}
				}
				// ���û�и�ָ�꣬��ֱ�Ӽ����Ԫ��
				if (!flag) {
					result.mat[i].add(new SCEntry(e2.col, new Complex(alpha * e2.val.real, alpha * e2.val.image)));
					result.nnz++;
				}
			}
		}
		return result;
	}
	
	/**
	 *     复矩阵加上一个实矩阵
	 * @param alpha
	 * @param B
	 * @return
	 */
	public SCMatrix add(double alpha, SRMatrix B) {
		SCMatrix result = new SCMatrix(this);
		int M = result.m;
		boolean flag = false;
		if (abs(alpha) == 0.0) {
			return result;
		}

		for (int i = 0; i < M; i++) {
			int mat2Size = B.mat[i].size();
			for (int j = 0; j < mat2Size; j++) {
				// flag��¼�Ƿ��ڵڶ��������д��ڸ�ָ��
				flag = false;
				SREntry e2 = B.mat[i].get(j);
				int resultSize = result.mat[i].size();
				for (int k = 0; k < resultSize; k++) {
					SCEntry e1 = result.mat[i].get(k);
					// Ѱ���������е���ͬ�����Ƿ�����ͬ����ָ��
					if (e1.col == e2.col) {
						double temp2Real = e1.val.real + alpha * e2.val;
						double t1 = Math.min(Math.abs(temp2Real), e1.val.norm());
						double ratio = t1 / sqrt(temp2Real * temp2Real + e1.val.image * e1.val.image);
						if (ratio > 1.0e-8) {
							e1.val.real = temp2Real;
						} else {
							result.mat[i].remove(k);
							result.nnz--;
						}
						flag = true;
						break;
					}
				}
				// ���û�и�ָ�꣬��ֱ�Ӽ����Ԫ��
				if (!flag) {
					result.mat[i].add(new SCEntry(e2.col, new Complex(alpha * e2.val, 0)));
					result.nnz++;
				}
			}
		}
		return result;
	}

	public void copy(SCMatrix A) {
		assert(m == A.m);
		assert(n == A.n);
		this.clear();
		for(int i = 0; i < m; i++) {
			for(SCEntry ele: mat[i]) {
				A.mat[i].add(new SCEntry(ele.col, ele.val));
				this.nnz++;
			}
		}
	}
	public void add(double alpha, SRMatrix B, SCMatrix result) {
		assert(m == B.getM());
		assert(n == B.getN());
		assert(m == result.m);
		assert(n == result.n);
		int M = result.m;
		boolean flag = false;
		result.copy(this);
		if (abs(alpha) == 0.0) {
			return;
		}

		for (int i = 0; i < M; i++) {
			int mat2Size = B.mat[i].size();
			for (int j = 0; j < mat2Size; j++) {
				// flag��¼�Ƿ��ڵڶ��������д��ڸ�ָ��
				flag = false;
				SREntry e2 = B.mat[i].get(j);
				int resultSize = result.mat[i].size();
				for (int k = 0; k < resultSize; k++) {
					SCEntry e1 = result.mat[i].get(k);
					// Ѱ���������е���ͬ�����Ƿ�����ͬ����ָ��
					if (e1.col == e2.col) {
						double temp2Real = e1.val.real + alpha * e2.val;
						if (sqrt(temp2Real * temp2Real + e1.val.image * e1.val.image) != 0.0) {
							e1.val.real = temp2Real;
						} else {
							result.mat[i].remove(k);
							result.nnz--;
						}
						flag = true;
						break;
					}
				}
				// ���û�и�ָ�꣬��ֱ�Ӽ����Ԫ��
				if (!flag) {
					result.mat[i].add(new SCEntry(e2.col, new Complex(alpha * e2.val, 0)));
					result.nnz++;
				}
			}
		}
	}

	public static SCMatrix Add(double alpha, SCMatrix A, double beta, SCMatrix B) {
		assert(A.m == B.m);
		assert(A.n == B.n);
		int m = A.m;
		int n = A.n;
		SCMatrix result = new SCMatrix(m, n);
		Complex[] temp = Complex.zeros(n);
		for(int i = 0; i < m; i++) {
			Complex.fill(temp, 0.0);
			for(SCEntry eleA: A.mat[i]) {
				temp[eleA.col].real += alpha * eleA.val.real;
				temp[eleA.col].image += alpha * eleA.val.image;
			}
			for(SCEntry eleB: B.mat[i]) {
				temp[eleB.col].real += beta * eleB.val.real;
				temp[eleB.col].image += beta * eleB.val.image;
			}
			for(int k = 0; k < n; k++) {
				if (temp[k].norm() != 0.0) {
					result.mat[i].add(new SCEntry(k, new Complex(temp[k].real, temp[k].image)));
					result.nnz++;
				}
			}
		}
		return result;
	}

	public static SCMatrix Add(double alpha, SCMatrix A, double beta, SRMatrix B) {
		assert(A.m == B.getM());
		assert(A.n == B.getN());
		int m = A.m;
		int n = A.n;
		SCMatrix result = new SCMatrix(m, n);
		Complex[] temp = Complex.zeros(n);
		for(int i = 0; i < m; i++) {
			Complex.fill(temp, 0.0);
			for(SCEntry eleA: A.mat[i]) {
				temp[eleA.col].real += alpha * eleA.val.real;
				temp[eleA.col].image += alpha * eleA.val.image;
			}
			for(SREntry eleB: B.mat[i]) {
				temp[eleB.col].real += beta * eleB.val;
			}
			for(int k = 0; k < n; k++) {
				if (temp[k].norm() != 0.0) {
					result.mat[i].add(new SCEntry(k, new Complex(temp[k].real, temp[k].image)));
					result.nnz++;
				}
			}
		}
		return result;
	}

	public static SCMatrix Add(double alpha, SRMatrix A, double beta, SCMatrix B) {
		assert(A.getM() == B.m);
		assert(A.getN() == B.n);
		int m = B.m;
		int n = B.n;
		SCMatrix result = new SCMatrix(m, n);
		for(int i = 0; i < m; i++) {
			for(SREntry eleA: A.mat[i]) {
				result.addElement(i, eleA.col, alpha * eleA.val);
			}
			for(SCEntry eleB: B.mat[i]) {
				result.addElement(i, eleB.col, new Complex(beta * eleB.val.real, beta * eleB.val.image));
			}
		}
		return result;
	}

	/**
	 * 
	 * @param d
	 * @return
	 */
	public SCMatrix mul(double d) {
		SCMatrix result = new SCMatrix(m, n);
		mmul(d, result);
		return result;
	}

	public void mul(double d, SCMatrix result) {
		ParamCheck.checkEqual(result.m, m);
		ParamCheck.checkEqual(result.n, n);
		mmul(d, result);
	}

	private void mmul(double d, SCMatrix result) {
		IntStream.range(0, m).forEach(i -> mat[i].forEach(ele -> {
			result.mat[i].add(new SCEntry(ele.col, ele.val.mul(d)));
			result.nnz++;
		}));
		/*
		 * for (int i = 0; i < m; i++) { for (SCEntry ele : mat[i]) {
		 * result.mat[i].add(new SCEntry(ele.col, ele.val * d)); result.nnz++; } }
		 */
	}
	
	public SCMatrix mul(Complex d) {
		SCMatrix result = new SCMatrix(m, n);
		mmul(d, result);
		return result;
	}

	public void mul(Complex d, SCMatrix result) {
		ParamCheck.checkEqual(result.m, m);
		ParamCheck.checkEqual(result.n, n);
		mmul(d, result);
	}

	private void mmul(Complex d, SCMatrix result) {
		IntStream.range(0, m).forEach(i -> mat[i].forEach(ele -> {
			result.mat[i].add(new SCEntry(ele.col, ele.val.mul(d)));
			result.nnz++;
		}));
		/*
		 * for (int i = 0; i < m; i++) { for (SCEntry ele : mat[i]) {
		 * result.mat[i].add(new SCEntry(ele.col, ele.val * d)); result.nnz++; } }
		 */
	}
	

	public Complex[] mul(Complex[] u) {
		ParamCheck.checkEqual(n, u.length);
		Complex[] v = Complex.zeros(m);
		mmul(u, v);
		return v;
	}

	public Complex[] mul(double[] u) {
		ParamCheck.checkEqual(n, u.length);
		Complex[] v = Complex.zeros(m);
		mmul(u, v);
		return v;
	}

	public void mul(Complex[] u, Complex[] v) {
		ParamCheck.checkEqual(n, u.length);
		ParamCheck.checkEqual(m, v.length);
		mmul(u, v);
	}

	private void mmul(Complex[] u, Complex[] v) {
		// Arrays.parallelSetAll(v, i -> mat[i].parallelStream().mapToDouble(ele ->
		// ele.val * u[ele.col]).sum());
		for (int i = 0; i < m; i++) {
			v[i].real = 0.0;
			v[i].image = 0.0;
			for (SCEntry ele : mat[i]) {
				Complex temp = ele.val.mul(u[ele.col]);
				v[i].real += temp.real;
				v[i].image += temp.image;
			}
		}
	}

	private void mmul(double[] u, Complex[] v) {
		// Arrays.parallelSetAll(v, i -> mat[i].parallelStream().mapToDouble(ele ->
		// ele.val * u[ele.col]).sum());
		for (int i = 0; i < m; i++) {
			v[i].real = 0.0;
			v[i].image = 0.0;
			for (SCEntry ele : mat[i]) {
				Complex temp = ele.val.mul(u[ele.col]);
				v[i].real += temp.real;
				v[i].image += temp.image;
			}
		}
	}

	public Complex[] mul(double coef, Complex[] u) {
		ParamCheck.checkEqual(n, u.length);
		Complex[] v = new Complex[m];
		mmul(coef, u, v);
		return v;
	}

	public void mul(double coef, Complex[] u, Complex[] v) {
		ParamCheck.checkEqual(n, u.length);
		ParamCheck.checkEqual(m, v.length);
		mmul(coef, u, v);
	}

	private void mmul(double coef, Complex[] u, Complex[] v) {
		// Arrays.parallelSetAll(v, i -> mat[i].parallelStream().mapToDouble(ele ->
		// ele.val * u[ele.col]).sum());
		for (int i = 0; i < m; i++) {
			v[i].real = 0.0;
			v[i].image = 0.0;
			for (SCEntry ele : mat[i]) {
				Complex temp = ele.val.mul(u[ele.col]);
				v[i].real += temp.real * coef;
				v[i].image += temp.image * coef;
			}
		}
	}

	public Complex[] mulT(Complex[] u) {
		ParamCheck.checkEqual(m, u.length);
		Complex[] v = new Complex[n];
		mmulT(u, v);
		return v;
	}

	public void mulT(Complex[] u, Complex[] v) {
		ParamCheck.checkEqual(m, u.length);
		ParamCheck.checkEqual(n, v.length);
		mmulT(u, v);
	}

	/**
	 *     矩阵的共轭转置乘以向量
	 * @param u
	 * @param v
	 */
	private void mmulT(Complex[] u, Complex[] v) {
		MVO.fill(v, 0.0);
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				Complex temp = ele.val.mulT(u[i]);
				v[ele.col].real += temp.real;
				v[ele.col].image += temp.image;
			}
		}
	}

	/**
	 * 
	 * @param rowDof
	 * @param colDof
	 * @param ele
	 */
	public void assemble(int[] rowDof, int[] colDof, Complex[][] ele) {
		int m = ele.length;
		int n = ele[0].length;
		ParamCheck.checkEqual(m, rowDof.length);
		ParamCheck.checkEqual(n, colDof.length);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				addElement(rowDof[i], colDof[j], ele[i][j]);
			}
		}

	}

	public void assemble(int[] rowDof, int[] colDof, double[][] ele) {
		int m = ele.length;
		int n = ele[0].length;
		ParamCheck.checkEqual(m, rowDof.length);
		ParamCheck.checkEqual(n, colDof.length);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				addElement(rowDof[i], colDof[j], ele[i][j]);
			}
		}

	}

	/**
	 * 
	 * @param dof
	 * @param ele
	 */
	public void assemble(int[] dof, Complex[][] ele) {
		assemble(dof, dof, ele);
	}
	public void assemble(int[] dof, double[][] ele) {
		assemble(dof, dof, ele);
	}

	public void assemble(int row, int[] colDof, Complex[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, colDof.length);
		for (int i = 0; i < m; i++) {
			addElement(row, colDof[i], ele[i]);
		}
	}

	public void assemble(int row, int[] colDof, double[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, colDof.length);
		for (int i = 0; i < m; i++) {
			addElement(row, colDof[i], ele[i]);
		}
	}

	/**
	 * 
	 * @param rowDof
	 * @param col
	 * @param ele
	 */
	public void assemble(int[] rowDof, int col, Complex[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, rowDof.length);
		for (int i = 0; i < m; i++) {
			addElement(rowDof[i], col, ele[i]);
		}
	}

	public void assemble(int[] rowDof, int col, double[] ele) {
		int m = ele.length;
		ParamCheck.checkEqual(m, rowDof.length);
		for (int i = 0; i < m; i++) {
			addElement(rowDof[i], col, ele[i]);
		}
	}



	/**
	 * 
	 */
	public void sort() {
		// for (int i = 0; i < m; i++) {
		// Collections.sort(mat[i], (a, b) -> a.col > b.col ? 1 : -1);
		// }
		IntStream.range(0, m).parallel().forEach(i -> Collections.sort(mat[i], (a, b) -> a.col > b.col ? 1 : -1));
		sorted = true;
	}

	/**
	 * 
	 * @return
	 */
	public boolean isSorted() {
		return sorted;
	}

	/**
	 * 
	 * @return
	 */
	public boolean makeEntrySymmetric() {
		int oldNNZ = nnz;
		for (int i = 0; i < m; i++) {
			int len = mat[i].size();
			for (int j = 0; j < len; j++) {
				int col = mat[i].get(j).col;
				if (getElement(col, i).norm() < Constants.Er) {
					setElementFast(col, i, new Complex(0.0, 0.0));
				}
			}
		}
		return oldNNZ != nnz;
	}

	public int[][] getRowRange() {
		int[][] rowRange = new int[m][2];
		sort();
		for(int i = 0; i < m; i++) {
			int size = mat[i].size();
			if (size > 0) {
				rowRange[i][0] = mat[i].get(0).col;
				rowRange[i][1] = mat[i].get(size - 1).col;
			} else {
				rowRange[i][0] = -1;
			}

		}
		return rowRange;
	}

	public int[][] getColRange() {
		int[][] colRange = new int[n][];
		for(int i = 0; i < n; i++) {
			colRange[i] = new int[]{m + 1, -1};
		}
		sort();
		for(int i = 0; i < m; i++) {
			for(SCEntry ele: mat[i]) {
				if(colRange[ele.col][0] > i) {
					colRange[ele.col][0] = i;
				}
				if(colRange[ele.col][1] < i) {
					colRange[ele.col][1] = i;
				}
			}
		}
		return colRange;
	}

	public SCMatrix mul(SRMatrix b) {
		assert(n == b.getM());
		int bn = b.getN();
		SCMatrix result = new SCMatrix(m, bn);
		for(int i = 0; i < m; i++) {
			for(SCEntry eleA: mat[i]) {
				for(SREntry eleB: b.mat[eleA.col]) {
					result.addElement(i, eleB.col, new Complex(eleA.val.real * eleB.val, eleA.val.image * eleB.val));
				}
			}
		}
		return result;
	}

	public String displayByCol() {
		StringBuilder sb = new StringBuilder();
		sb.append(m + "\t" + n + "\t" + nnz + "\t" + false + "\n");
		ArrayList<Pair<Integer, Complex>>[] cscmat = new ArrayList[n];
		for(int i = 0; i < n; i++) {
			cscmat[i] = new ArrayList<>();
		}

		for(int i = 0; i < m; i++) {
			for(SCEntry ele: mat[i]) {
				cscmat[ele.col].add(new Pair<>(i, ele.val));
			}
		}
		IntStream.range(0, n).parallel().forEach(i -> Collections.sort(cscmat[i], (a, b) -> a.first > b.first ? 1 : -1));
		for(int i = 0; i < n; i++) {
			for(Pair<Integer, Complex> ele: cscmat[i]) {
				sb.append(ele.first + "\t  " + i + "\t  " + ele.second + "\n");
			}
		}
		return sb.toString();
	}
	@Override
	public String toString() {
		StringBuilder sb = new StringBuilder();
		sb.append(m + "\t  " + n + "\t  " + nnz + "\t  " + isHermitian() + "\n");
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				sb.append((i + 1) + "\t  " + (ele.col + 1) + "\t  " + ele.val + "\n");
			}

		}
		return sb.toString();
	}

	public String display() {
		StringBuilder sb = new StringBuilder();
		sb.append(m + "\t" + n + "\t" + nnz + "\t" + isHermitian() + "\n");
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				sb.append(i + "\t" + ele.col + "\t" + ele.val + "\n");
			}

		}
		return sb.toString();
	}

	public Complex[] diag() {
		Complex[] d = new Complex[m];
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				if (ele.col == i) {
					d[i].real = ele.val.real;
					d[i].real = ele.val.image;
				}
			}
		}
		return d;
	}

	public SCMatrix transpose() {
		SCMatrix result = new SCMatrix(n, m);
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				result.mat[ele.col].add(new SCEntry(i, ele.val));
			}
		}
		return null;
	}

	public SCMatrix triU() {
		SCMatrix upper = new SCMatrix(m, n);
		for (int i = 0; i < m; i++) {
			int ii = i;
			mat[i].stream().filter(ele -> ele.col > ii).forEach(ele -> upper.setElementFast(ii, ele.col, ele.val));
		}
		return upper;
	}

	public SCMatrix triL() {
		SCMatrix lower = new SCMatrix(m, n);
		for (int i = 0; i < m; i++) {
			int ii = i;
			mat[i].stream().filter(ele -> ele.col < ii).forEach(ele -> lower.setElementFast(ii, ele.col, ele.val));
			// lower.mat[i].addAll(mat[i].stream().filter(ele -> ele.col < ii));
		}
		return lower;
	}

	/**
	 * ʹ����Ԫ�صõ������Ԫ�ص�ֵ
	 * 
	 * @param ii��������ָ��
	 * @param jj��������ָ��
	 * @param values������Ԫ��ֵ
	 */
	public void getTriValues(int[] ii, int[] jj, Complex[] values) {
		int in = 0;
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				ii[in] = i;
				jj[in] = ele.col;
				values[in] = new Complex(ele.val);
				in++;
			}
		}
	}

	/**
	 * ����õ�ϡ�����Ԫ�����ֵ @return��ϡ�����Ԫ�����ֵ
	 */
	public double maxNormElement() {
		double[] maxval = new double[m];
		for (int i = 0; i < m; i++) {
			maxval[i] = mat[i].stream().mapToDouble(ele -> ele.val.norm()).max().orElse(0.0);
		}
		return Arrays.stream(maxval).max().orElse(0.0);
	}

	/**
	 * ����õ�ϡ�����Ԫ�����ֵ�����������ֵ��������ֵ
	 * 
	 * @param IJ�����ֵ��������ֵ @return��ϡ�����Ԫ�����ֵ
	 */
	public Complex maxNormElement(int[] IJ) {
		Complex[] maxval = new Complex[m];
		int[] maxJ = new int[m];
		for (int i = 0; i < m; i++) {
			if (mat[i] != null) {
				int tempJ = mat[i].get(0).col;
				Complex tempVal = mat[i].get(0).val;
				for (int j = 1, size = mat[i].size(); j < size; j++) {
					SCEntry ele = mat[i].get(j);
					if (ele.val.norm() > tempVal.norm()) {
						tempJ = ele.col;
						tempVal = ele.val;
					}
				}
				maxJ[i] = tempJ;
				maxval[i] = tempVal;
			}
		}
		Complex result = maxval[0];
		IJ[0] = 0;
		IJ[1] = maxJ[0];
		for (int i = 1; i < m; i++) {
			if (result.norm() < maxval[i].norm()) {
				result = maxval[i];
				IJ[0] = i;
				IJ[1] = maxJ[i];
			}
		}
		return result;
	}

	public boolean isHermitian() {
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[i]) {
				Complex temp = getElement(ele.col, i);
				if (abs(ele.val.real - temp.real) > Constants.Er || abs(ele.val.image + temp.image) > Constants.Er) {
					return false;
				}
			}
		}
		return true;
	}
	
	

	public int[] symamd() {
//			if(isSymmetric()) {
//				int[] degree = new int[m];
//				Arrays.setAll(degree, i -> mat[i].size());
//				System.out.println(Arrays.toString(degree));
//				int[] perm = new int[m];
//				Algorithms.selectionSort(degree, perm, true);
//				System.out.println(Arrays.toString(perm));
//			} else {
//				System.out.println("The matrix is not symmetric");
//	//			return null;
//			}
//		int[] degree = new int[m];
//		boolean[] flag = new boolean[m];
//		Arrays.setAll(degree, i -> mat[i].size());
//		Pair<Integer, Integer> info = new Pair<>();
//		ArrayList<Integer> idxList = new ArrayList<>();
//		for(int i = 0; i < m; i++) {
//			if (flag[i] == false) {
//				for(SCEntry ele : mat[i]) {
//					if (flag[ele.col] == false) {
//						
//					}
//				}
//			}
//			//for(int j = i; j < m; j++) {
//			//	
//			//}
//		}
//		
//		Arrays.setAll(degree, i -> mat[i].size());
//		Pair<Integer, Integer> info = new Pair<>();
//		
//		MVO.min(degree, info);
//		int idx = info.getFirst();
//		for(SCEntry ele: mat[idx]) {
//			flag[ele.col] = true;
//		}

		int[] degree = new int[m];
		int[] perm = new int[m];
		boolean[] flag = new boolean[m];
		Arrays.setAll(degree, i -> mat[i].size());
		Pair<Integer, Integer> info = new Pair<>();
		for (int i = 0; i < m; i++) {
			MVO.min(degree, info);
			int idx = perm[i] = info.getFirst();
			for (SCEntry ele : mat[idx]) {
				if (flag[ele.col] == false) {
					degree[ele.col]--;
				}
			}
			flag[idx] = true;
			degree[idx] = Integer.MAX_VALUE;
		}
		/*
		 * System.out.println(Arrays.toString(perm)); flag = new boolean[m]; for(int i =
		 * 0; i < m; i++) { if (flag[i] == false && i != perm[i]) {
		 * System.out.println("change: " + i + "\t, " + perm[i]); this.exchangeRow(i,
		 * perm[i]); System.out.println(MVO.toString(this.toDenseMatrix()));
		 * this.exchangeCol(i, perm[i], nBand);
		 * System.out.println(MVO.toString(this.toDenseMatrix()));
		 * System.out.println("=================="); flag[i] = true; flag[perm[i]] =
		 * true; } }
		 */

		// int[] inversePerm = new int[m];
		// Arrays.setAll(inversePerm, i -> perm[i]);

		// System.out.println(Arrays.toString(perm));
		SCMatrix temp = new SCMatrix(m);
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : mat[perm[i]]) {
				temp.mat[i].add(new SCEntry(ele.col, ele.val));
			}
		}

		// System.out.println(MVO.toString(temp.toDenseMatrix()));
		int[] inversePerm = new int[m];
		for (int i = 0; i < m; i++) {
			inversePerm[perm[i]] = i;
		}
		// System.out.println(Arrays.toString(inversePerm));
		this.clear();
		for (int i = 0; i < m; i++) {
			for (SCEntry ele : temp.mat[i]) {
				mat[i].add(new SCEntry(inversePerm[ele.col], ele.val));
			}
		}
		temp.clear();
		return perm;
	}

	public DZcs_common.DZcs toDZcs_matrix() {
		DZcs_common.DZcs T;
		T = DZcs_util.cs_spalloc(0, 0, 1, true, true); /* allocate result */
		String line;
		for(int k = 0; k < m; k++) {
			for(SCEntry ele: mat[k]) {
				if (!DZcs_entry.cs_entry(T, k, ele.col, new double[]{ele.val.real, ele.val.image}))
					return (null);
			}

		}
		return (T);
	}
}
